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PREFACE 



Several years ago it was suggested that I should deliver 
two lectures on the Monte Carlo method to the students of 
the Computer Technology Faculty of the Public University. 
I agreed. On the eve of the first lecture I found out to my 
dismay that most of the audience were not familiar with 
the probability theory. It was too late to retreat. I had no 
other choice but to present in the course of the lecture a com- 
plementary section which surveyed the fundamentals of the 
probability theory. 

These lectures were repeated several years running. Gra- 
dually their content has "set". The complimentary para- 
graph (Sec. 2) is retained in this edition, and I feel that it 
merits a few words on my part. 

Everybody had at some moment used the words "proba- 
bility" and "random variable". The intuitive idea of the 
probability (considered as frequency) corresponds more or 
less to the true meaning of this concept. But as a rule the 
intuitive idea of the random variable differs quite conside- 
rably from the mathematical definition. Thus, the notion 
of the probability is assumed known in Sec. 2, and only the 
more complicated notion of the random variable is clari- 
fied. This section cannot replace a course in the probability 
theory: the presentation is simplified and proofs are omitted. 
But it still presents certain concept of the random variables 
sufficient for understanding of Monte Carlo techniques. 
The basic aim of this book is to prompt the specialists 
in various branches of knowledge to the fact that there are 
problems in their fields that can be solved by the Monte 
Carlo method. 

The problems considered in these lectures are sufficiently 
simple and diverse. Of course, they could not encompass all 
the problems this method may be applied to. One example 
will be sufficient to illustrate this statement. Not a word 



in this book is devoted to medicine. But the methods de- 
scribed in Sec. 7 make it possible to calculate irradiation doses 
in radio-therapy. If we have a program for computing radia- 
tion absorbed in various body tissues we can select both 
the dose and direction of irradiation most efficiently and 
with no damage done to normal tissues. 

Everything that was presented in the lectures is included 
into this text, the examples being given in some more 
detail. An additional section (Sec. 9) is provided. 

/. Sobol 
Moscow, 1967 

The second edition is practically identical to the first. 
Only negligible editors' corrections were introduced. 

/. Sobol 
Moscow, 1971 



Introduction 



See. 1. General 

The Monte Carlo method is a numerical method of solving 
mathematical problems by means of random sampling. 

1.1. The Origin of the Monte Carlo Method. The method 
is conventionally assumed to be born in 1949 when the paper 
"The Monte Carlo method" was published*. American mathe- 
maticians John von Neumann and Stanislav Ulam are thought 
to have developed this method. First papers on the Monte 
Carlo method in the USSR were published in 1955-1956**. 

Incidentally, the theoretical foundation of the method 
has been known already for a long time. Furthermore, cer- 
tain problems in statistics were sometimes calculated by 
means of random sampling, i.e. in fact by means of the 
Monte Carlo method. Until the advent of electronic compu- 
ters, however, this method could not be used on any signi- 
ficant scale because manual simulation of random variables 
is a very labour-consuming procedure. Thus, the Monte 
Carlo method as an extremely versatile numerical technique 
became feasible only owing to the computers. 

As for the term "Monte Carlo", it is due to the name of 
the city of Monte Carlo in the Principality of Monaco famous 
for its casino. The point is that the roulette is one of the 
simplest mechanical devices for ... generation of random 
numbers. We will discuss it in Sec. 3. But it appears worth- 
while to answer here one frequently posed question: "Does 
the Monte Carlo method help to win in roulette?" The ans- 
wer is "no". It is not even concerned with it. 



* Metropolis N., Ulam S. The Monte Carlo method, /. Amer. 
statistical assoc, 1949, 44, No. 247, 335-341. 

** Those were the papers by V. V. Chavchanidze, Yu.A. Shreider 
and V.S. Vladimirov. 



1.2. Example. To help the reader get a clear idea of what 
it is all about we will take a fairly simple example. Let us 
assume that we have to calculate the surface area of a plane 
figure S. It may be a completely arbitrary figure with 
curvilinear boundary whether it be connected or consisting 




FIG. 1 

of several sections, and specified graphically or analytically. 
Let it be a figure drawn in Fig. 1 and let us assume that it 
is enclosed within a unit square. 

Let us select N random points inside the square. Let N' 
designate the number of points that happened to fall within 
S. It is geometrically obvious that the area S is appoximate- 
ly equal to the ratio N'IN. The greater N the higher the 
accuracy of the estimate. 

The number of points selected in Fig. 1 is TV = 40. Of 
these points N' = 12 are inside S. The ratio N'IN = 
= 12/40 = 0.30 while the true area S is 0.35*. 



* In practice the Monte Carlo method is not applied to calculate 
the area of the plane figure because other methods are available which, 
though more complicated, provide much higher accuracy. 

However, the Monte Carlo method, as illustrated in the example 
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1.3. Two Distinctive Features of the Monte Carlo Methud. 

The first feature is the simple structure of the computation 
algorithm. As a rule, a program is prepared to perform only 
one random trial (in Example 1.2 we have to select a random 
point within the limits of the square and check whether it 
lies within S). This trial is repeated N times, each trial 
being independent of all others, and the results of all trials 
are averaged. 

Thus, the Monte Carlo method is sometimes called the 
method of statistical trials. 

The second feature of the method is as follows: the error 
of calculations is, as a rule, proportional to Y DIN where D 
is a constant and N is the number of trials. This formula shows 
that 100-fold increase of N (and thus of the work) is neces- 
sary to reduce the error by a factor of 10 (in other words, 
to obtain one more accurate decimal digit in the result). 

It is clear that high accuracy cannot be attained with 
such an approach. Thus it is usually said that the Monte 
Carlo method is especially efficient in solving those problems 
which do not require a high degree of accuracy (e.g. 5 to 
10 per cent). 

However, a given problem can be solved through various 
versions of the Monte Carlo method*, each having a diffe- 
rent value of D. Accuracy can be considerably improved 
for many problems by a proper choice of the computation 
technique with an appreciably smaller value of D. 

1.4. Problems Solved by the Monte Carlo Method. First 
of all, the Monte Carlo method allows to simulate any pro- 
cess whose development is affected by random factors. Se- 
condly, many mathematical problems not affected by any 
random influences can be connected with an artificially 
constructed probabilistic model (even more than one) mak- 
ing possible the solution of these problems. This in fact 
has been done in Example 1.2. 



above, permits to calculate just as simply the "multidimensional 
volume" of a body in a multidimensional space. In this case the Monte 
Carlo method is often the only numerical method enabling us to solve 
the problem. 

* It becomes more widespread nowadays in publications abroad 
to speak of Monte Carlo methods (in plural), meaning that one and 
the same problem can be computed by means of simulating various 
random variables. 
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Thus wo can consider the Monte Carlo technique as a 
universal method of solving mathematical problems. 

Of particular interest is the fact that in certain cases it 
is expedient to use an artificial model instead of simulat- 
ing a real random process. Such a problem is discussed in 
Sec. 7. 

1.5. More about the Example. Let us go back to Example 
1.2. The computation required the selection of random points 
within the unit square. How should we actually realize it? 




FIG. 2 

Let us imagine the following experiment. Fig. 1 scaled 
up with the figure S and the square is fixed on the wall as a 
target. A marksman, placed some distance from the wall, 
shoots N times aiming at the centre of the square. Of course, 
the bullets will not all strike the centre of the target: they 
will pierce the target in N random points*. Can we estimate 
area S from these points? 



* We assume that the marksman is not the world champion and is 
removed at a sufficiently long distance from the target. 
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The result of suck an experiment is shown in Fig. 2. In 
this experiment N = 40, N' = 24 and the ratio iV'AV = 
= 0.60 which is almost twice the true value of the area 
(0.35). It is quite obvious however that very high quali- 
fication of the marksman will lead to a very bad result of 
the experiment, because almost all punctures will lie close 
to the centre, i.e. inside S*. 

It is easily seen that the above method of calculating 
the area will be valid only if the random points are not 
"just random" but also "uniformly scattered" over the whole 
square. In order to attach a precise meaning to these terms 
we have to introduce the definition of random variables and 
to cite some of their properties. This information is given 
without proof in Sec. 2. Readers familiar with any course 
in the probability theory can skip the whole of Sec. 2 with 
the exception of Item 2.5. 



* The manner in which random points were selected in Figs. 1 
and 2 will be elucidated in Sec. 4.5. 



CHAPTER I 

Simulation of Random 
Variables 

Sec. 2. Random Variables 

We shall assume that the reader is more or less familiar 
with the concept of probability and shall turn directly to 
the concept of random variable. 

The term "random variahle", in its commonplace sense 
is used to emphasize that one does not know what specific 
value this variable will assume. Sometimes these words 
merely hide the fact that one simply docs not know the value 
of this variable. 

Mathematicians, on the contrary, use the same expression 
"random variahle" in quite a definite sense. They say that 
though we indeed do not know what value will this variable 
take on in a given case, we do know what values it can as- 
sume and what are the respective probabilities. The result 
of a single trial representing this random variahle cannot 
be accurately predicted from these data but we can predict 
very reliably the result of a large number of trials. The grea- 
ter the number of trials (or, as it is normally expressed, 
the larger the statistics) the more accurate will our predic- 
tions be. 

Thus, to present a random variable we must indicate 
what values it can assume and what are the probabilities 
of these values. 

2.1. Discrete Random Variables. If a random variable £ 
takes on values from a discrete set x 1 , x 2 , . . ., x n it is 
called a discrete random variable.* 

The discrete random variable | is specified by the table 

(X i x<> . . . x n \ 
ftft-J- (T > 

* The probability theory also considers discrete random variables 
which can assume infinite number of values. 
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where x lt a\>, . . ., x n are tlie possible values of the variable 
E and Pi, p 2 , . . ., p n are the corresponding probabilities. 
To be precise, the probability that the random variable | 
will be equal to x t (we will designale it as P {E = xi)) is 
equal to p t : 

P {E = x,} = Pi . 

The table (T) is called the distribution of the random va- 
riable. 

Generally speaking, the values x l7 x 2 , . ■ ., x n can be 
arbitrary. But the probabilities p x , p 2 , . . ., p n must 
satisfy two conditions: 

(a) all Pi are positive: 

Pi > 0; (1) 

(b) the sum of all p t is equal to 1: 

Pi + Pa + • ■ • + Pn = I- (2) 

The last condition means that in each trial E necessarily 
takes on one of the values x x , x 2 , . . ., x n . 
The quantity 

Mt-£x iPi (3) 

is called the mathematical expectation, or the expected value 
of the random variable E. 

To elucidate the physical meaning of this value we will 

rewrite il in the following form: 



ME- 



As follows from this relation, ME is the mean value of the 
variable E and the values of x t associated with greater pro- 
babilities are included into the sum with larger weights*. 



* Weight-averaging is encountered very often and in most different 
spheres of science. In mechanics, for example, if masses m u m 2 , . . ., m n 
are located at the points x lt x 2 , . . ., x n (on the Ox-axis), then the 
abscissa of the centre of gravity of this system is calculated by means 
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H 

~2jPi 



Let us discuss the basic properties of the expected value. 
If c is an arbitrary non-random number, then 

M (| + c) = Mg + c, (4) 

M (c|) = cM£; (5) 

if i and n. are two arbitrary random variables, then 

M (i + n) = ME + Mr). (6) 
The quantity 

Di = M 1(1 - Mm (7) 

is called the variance of the random variable |. Thus the 
variance Di is the expected value of the square of the devia- 
tion of the random variable £ from its mean value M|. 
It is obvious that in all cases Di > 0. 

Expected value and variance are the two most important 
numerical characteristics of the random variable |. What 
is their practical significance? 

If in the course of many observations of the variable i 
we obtain a set of values i 1( i 2 , . . ., i N (each of these 
will be equal to one of the numbers x ly x. 2 , . . ., x n ), then 
the arithmetic mean of these numbers will be close to Mi 

-j^Gi +£•+■•■+!*)« Ml. (8) 

The variance D| characterizes the spreading of these quan- 
tities around the mean value Mi. 

Formula (7) for the variance can be transformed by means 
of Eqs. (4)-(6): 

Di = MU 2 -2Mi-i + (Mi) 2 l = 

= M(i 2 ) -2Mi-Mi + (Mi) 2 , 



whence 



Di = M (| 2 ) - (MH) 2 . (9) 



of the following formula: 



2j x i m i 
i=\ 
x = - 



n 

Obviously, in this case the sum of masses should not necessarily 
equal unity. 
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Usually the calculation of variance via Eq. (9) is simpler 
than by Eq. (7). 

Let us mention the basic properties of variance: if c 
is an arbitrary non-random number, then 

d (i -i- c) - m, (io) 

D (eg) = cDE. (11) 

The concept of independence of random variables plays 
an important role in the probability theory. This is in fact 
a rather complicated concept, though it may be quite clear 
in the simplest cases*. 

The following relations are valid for independent random 
variables E and T) 

M(|ti) = ME -Mi], (12) 

D (I + n) = D| + Dt|. (13) 

Example. Let us consider a random variable x with the distri- 
bution specified by the table 



1 


2 


3 


4 


5 


(i 


1 


1 


1 


1 


1 


1 


6 


G 


6 


6 


(i 


6 



It is obvious that the number of points obtained in one throw of a die 
can be considered a realization of this variable: any of the values is 
equally probable. Let us calculate the expected value and the variance 
of x. According to (3) 

Mx=l-4+--- -•f-6~=3.5. 
6 (> 

According to (9) 

Dx = M(x 2 ) — (Mx) 2 = l=. •! + ... +6 2 -i— (3.5) 2 = 2.917. 

Example. Let us consider a random variable 6 with the distri- 
bution 

/3 4 

B= 1 1 

\ 2 2 

The game of tossing a coin, or "heads or tails", with the condition 
that a toss of heads brings 3 points and that of tails, 4 points, can be 

* Let us assume simultaneous observation of random variables | 
and T). If the distribution of the variable § does not change after we 
have found the value assumed by the variable n, then it is natural to 
consider % independent of n. 
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considered a realization of this variable. In this case 
M(6) = 0.5-3 + 0.5-4 = 3.5; 
D (6) = 0.5 (32 + 42) - (3.5) 2 = 0.25. 

We see that M9 = Mx, but DO < Dx. The inequality could be 
easily predicted because the maximum value of deviation of from 
3.5 is +0.5 while the spreading of x can reach +2.5. 

2.2. Continuous Random Variables. Let us assume that 
a certain amount of radium is placed on a plane at the origin 
of coordinates. An a-particle is emitted in each decay of 
a radium atom. The direction of motion of this particle 




FIG. 3 



will be characterized by an angle \p (Fig. 3). Since^theoreti- 
cally and practically any direction^ of emission is possible, 
this random variable can assume any value from 
to 2k. 

We shall call a random variable % continuous if it takes 
on any value out of a certain interval (a, b). 

A continuous random variable | is defined by specifying 
the interval (a, b) of its variation, and the function p (x), 
called the probability density of the random variable | 
(or the distribution density of £). 

The physical meaning of p (x) is the following: let (a', b') 
be an arbitrary interval contained within (a, b) (i.e. a ^ 
^ a', b' sg; b). Then the probability that | falls inside 
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(a 1 , b') is given by the integral 



'{«'<!<&'}= j p(x)dx. 



(14) 



This integral (14) is equal to the hatched area in Fig. 4. 

The set of values of £ can be any interval. The cases when 
a = — oo as well as b — oo are also possible. However the 
density p (x) must satisfy two conditions similar to condi- 
tions (1) and (2) for discrete variables: 

(a) the density p (x) is positive: 

P (x) > 0; (15) 




(b) the integral of the density p (x) over the whole inter- 
val (a, b) is equal to 1: 



The value 



u 

\ p (x) dx — - 

a 

b 

l|= jxp(x) 



1. 



dx 



(16) 



(17) 



is called the expected value of the continuous random vari- 
able. 

This characteristic has the same meaning as in the case 
of a discrete random variable. In fact, it directly follows 
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from 

b 

\ xp (x) dx 

I p (x) dx 
a 

that this is the mean value of £: indeed, any value of x from 
the interval (a, b), entering the integral with the weight 
p (x) can be the value of |*. 

All the results of Sec. 2.1, starting with Eq. (4) and inclu- 
ding Eq. (13) are valid for continuous random variables. 
The definition of variance (7) and the computation formula 
(9) as well as all the statements with respect to properties 
of Mi and D£ are also valid. We shall not repeat them. 

We will mention only one more expression for the expec- 
ted value of a function of %. Let us assume again that a ran- 
dom variable £ is characterized by the probability density 
p (x). Let us take an arbitrary continuous function / (x) 
and introduce a random variable n = / (|). It can be shown 
that 

6 

M/(|)= y!(x)p(x)dx. (18) 

a 

We must ^emphasize that in general M/ (£) =£ / (M£). 

A random variable y, defined on the interval (0, 1) with 
the density p (x) = 1 is said to be uniformly distributed 
in (0, 1) (Fig. 5). 

Indeed, for any interval (a', b') within (0, 1) the proba- 
bility that y will take on the value within (a', b') is equal to 

b' 

\ p (x) dx — b' — a', 



* In this case it is equally possible to indicate a similar formula 
from mechanics: if the linear density of a rod a ^ x ^ b is p (x), then 
the abscissa of the rod gravity centre is calculated by means of the 
expression 

b 

\ xp (x) dx 



b 

J p (x) dx 
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i.e. to the length of this interval. If, in particular, we divide 
the (0, 1) interval into an arbitrary number of intervals 
of equal length, the probability of 7 falling within any one 
of these intervals will be the same. 



Yl\ 



P(*) 



FIG. 5 

It can be easily shown that 

1 1 



My = I xp(x)dx — \ x dx --- -5- ; 


1 

Dy = \ x 2 p (x) dx — (My) 2 = -J - T 
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This random variable y will be more than once encountered 
in further discussion. 

2.3. Normal Random Variables. A random variable £» 
defined on the entire axis ( — 00, 00) and characterized by 
the density 

(«~q)8 

p(x)--=- 7 L—e' ™ , (19) 

where a and a > are numerical parameters*, is said to be 
normal (or Gaussian) random variable. 

The parameter a does not affect the shape of the p (x) 
curve; its variation only shifts the curve as a whole along 

* a is a number and not a random variable. The Greek letter is 
used here as the generally accepted one. 
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the z-axis. On the contrary, variation of a changes the shape 
of the curve. It is indeed easily seen that 

max p (x) - p (a) = J . 
a y in 

A decrease of a is accompanied by an increase of max p (x). 
But according to the condition (16) the entire area below 
the p (x) curve is equal to unity. Therefore the curve will 
stretch upward in the vicinity of x = a but will decrease 



(5=0.5 




FIG. 6 

at all sufficiently large values of x. Two normal densities 
corresponding to a = 0, a — 1 and a — 0, a = 0.5 are 
plotted in Fig. 6. (Another normal density is plotted in 
Fig. 21, see below, p. 49.) 
It can be proved that 

MS - a, DC = 2 . 

Normal random variables are frequently encountered in 
investigation of most diverse problems. The reason will 
be discussed somewhat later. For example, the experimental 
error 8 is usually a normal random variable. If there is no 
systematic error (bias) then a — M8 = 0. The value o — 
= \^D8, called the standard deviation, characterizes the 
error of the method of measurement. 
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The "three sigma" rule. It is easily calculated that regard- 
less of the values of a and a in (19) 

a+Sa 

\ p(x)dx = 0.997. 

a- 3a 

As follows from (14), 

P {a — 3a < £ < a + 3a} = 0.997. (20) 

The probability 0.997 is so close to 1 that sometimes the 
last formula is given the following interpretation: it is 
practically impossible to obtain in a single trial a value of Z, 
deviating from M£ by more than 3a. 

2.4. Central Limit Theorem of the Probability Theory. 
The statement of this spectacular theorem has first been 
given by Laplace. A number of outstanding mathematicians, 
P. L. Chebyshev, A. A. Markov and A. M. Lyapunov among 
them, investigated the problems of generalizing this theo- 
rem. Its proof is rather complicated. 

Let us consider TV identical independent random variab- 
les \ x , i 2 , . . ., Ijv such that their probability distribu- 
tions coincide. Consequently, both their expected values and 
variances coincide as well. 

Let us designate 

Mgj = M£ 2 = . . . = MIjv = m, 
D^ = D£ 2 = . . = Dl„ = b\ 

and designate the sum of all these values as p N 

Pn = li + l 2 + • ■ • + In- 

It follows from Eqs. (6) and (13) that 

Mp w = iA(li + U + ■ ■ ■ + In) = Nm, 
Dp N = Dfo + I, + . . . + l N ) = Nb\ 

Let us consider now a normal random variable t, N with 
the same parameters: a = Nm, a 2 = Nb 2 . 

The central limit theorem states that for any interval 
{a', b') for large N 

b' 

P{a'<p w <6'}» j Pt N (z)dx. 

a' 
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The physical meaning of this theorem is obvious: the 
sum p N of a large numher of identical random variables is 
approximately normal (p PN (x) « p^ N (*))■ 

Actually this theorem is valid for much more general 
conditions: all the addends £ x , £ 2 , . . ., | A - are not requi- 
red to be identical and independent, it is only essential 
that none of them play too great a role in this sum. 

It is this theorem that explains why normal random varia- 
bles are met so often in nature. And indeed, each time we 
come across an aggregate effect of a large number of negli- 
gible random factors we find that the resulting random 
variable is normal. 

Deviation of an artillery shell off the target, for example, 
is almost always a normal random variable because it 
depends on meteorological conditions in various sections 
of the trajectory as well as on many other factors. 

2.5. The General Scheme of the Monte Carlo Method. Let 
us assume that we need to calculate some unknown value m. 
We shall try to find a random variable % such that M£ = m. 
Let us assume therewith that D| = b 2 . 

Let us consider N random independent variables t ly | 2 , . . . 
. . ., £jv with distributions identical to that of H. If N is 
sufficiently large then it follows from 2.4 that the distribu- 
tion of the sum p N =- | x + £ 2 + • • • + In will 1> C appro- 
ximately normal and will have the parameters a = Nm, 
o 2 = Nb*. Formula (20) yields 

F{Nm - 36 VW< p N < Nm + 3b V N} « 0.997. 

If we divide the inequality in braces by N we will arrive 
at an equivalent inequality and the probability will not 
change: 

P / m » <Pw <m + _|L\« 0.997. 

We will rewrite the last expression in a somewhat different 
form: 

p{|-fS^H<wH - 997 - (21) 



;'=i 



This formula is very important to the Monte Carlo method. 
It gives us both the method of calculating m and the esti- 
mate of the error. 
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Let us indeed sample N values of the random variable |*. 
The formula (21) indicates that the arithmetic mean of 
these values will be approximately equal to m. The error 
of this approximation will most probably not exceed the 
value 3b/y~N. This error evidently approaches zero as 
N increases. 

Sec. 3. Generation of Random Variables 
by Electronic Computers 

Sometimes the statement "generation of random varia- 
bles" by computers is the cause of bewilderment: everything 
performed by a computer must be programmed in advance, 
so where does the chance come in? 

This statement is indeed fraught with certain difficulties 
but these mostly concern philosophy and we shall not con- 
sider them here. 

We will only note that random variables dealt with in 
Sec. 2 represent an ideal mathematical concept. Whether 
some natural phenomenon can be actually described by 
means of these variables can only be found out experimen- 
tally. Such a description is always approximate. Moreover, 
a random variable which describes quite satisfactorily 
a certain physical quantity in one type of phenomena may 
prove quite an unsatisfactory characteristic of the same 
quantity when other phenomena are being investigated. 

The situation is similar to that of a road which can be 
considered straight on the map of the country (ideal mathe- 
matical straight line with "zero width"), but appears as a 
twisted band on a large-scale map of a settlement... 

The methods of generating random variables are usually 
divided into three types: tables of random numbers, 
random number generators and the method of pseudo-ran- 
dom numbers. 

3.1. Random Number Tables. Let us run the following 
experiment. Let us write down the digits 0, 1, 2, . . ., 9 
on ten identical sheets of paper, put them into a hat, mix 
them and then extract one sheet at a time, each time return- 
ing it to 1he hat and again mixing the sheets. We shall list 

* Determining a single value of each of the variables | l , £,_> • • •< Ijv 
is equivalent to determining N values of a single variable §, because 
all these random variables are identical (have identical distributions). 
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the numbers obtained in this way in a table similar to Table 
A at the end of this book on page 74. (For the sake of conve- 
nience the numbers in Table A are arranged in groups of 
5 digits.) 

Such a table is called the Random Number Table, though 
"The Table of Random Digits" would be a better title. This 
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table can be introduced into the memory of a computer. 
When later in the course of calculations we shall need values 
of the random variable with the distribution 




0.1 



1 

0.1 



2 
0.1 



9 
0.1 



(22) 



we shall retrieve consecutive numbers from this table. 
The most comprehensive of all published tables of ran- 
dom numbers contains 1000000 numbers*. Of course, this 
table has been compiled by means of a device more advanced 
than a hat: a special roulette with electric devices was 
constructed. The simplest design of such a roulette is shown 
in Fig. 7 (a rotating disk is abruptly stopped and the number 
indicated by a fixed pointer is recorded). 

* RAND Corporation, A Million Random Digits with 100 000 
Normal Deviates, The Free Press, 1955. 
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It should be noted that a good table of random numbers 
is not so easy to compile as one might expect. Any real 
physical instrument generates random numbers with distri- 
bution that to a certain extent differs from the ideal distri- 
bution (22). Besides, errors are possible during the experi- 
ment (for instance, one of the paper sheets may stick for 
some time to the lining of the hat). For this reason the com- 
piled tables are thoroughly checked by means of special 
statistical tests: whether some of the properties of the group 
of numbers contradict the hypothesis that these numbers 
are the values of the random variable (22). 

Consider one of the simplest tests. Let us take a table 
containing N numbers, in which v is the number of zeros, 
v t is the number of unities, v 2 — that of two's, and so on. 
Let us calculate the sum 

S(v,- 0.17V) 2 . 

i=0 

The probability theory enables us to predict in what range 
this sum should be found: its value should not be too large 
(since the expected value of each v ; is equal to 0.1A) but 
neither should it be too small (because it would mean "too 
regular" distribution of numbers). 

Random number tables are used in Monte Carlo calcula- 
tions only if these are performed manually. The thing is, 
all computers have a comparatively small internal store, 
and a large table cannot be stored in it; storage of the table 
in the external data store and constantly retrieving them 
considerably slows down the computation. 

It is not excluded though that time will bring a sharp 
increase of memory capacity in computers; this will make 
the use of random number tables quite practicable. 

3.2. Generators of Random Numbers. It might seem that 
a roulette mentioned in Sec. 3.1 could be coupled to a com- 
puter and would generate random numbers as soon as they 
are required. However, any mechanical device will be too 
slow for a computer. Therefore, noise signals in electron 
tubes are most often used as generators of random varia- 
bles: if the noise level has exceeded a fixed threshold an even 
number of times during a certain fixed time interval At, zero 
is recorded, and if this number is odd, unity is recorded.* 

* There also exist more advanced devices. 
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At a first glance this appears to be a very convenient 
method. Let m such generators operate simultaneously and 
send random zeros and unities into all binary digit places 
of a special cell. Each cycle produces one m-digit number. 
At any phase of computations we can refer to this cell and 
retrieve a value of the random variable y distributed uni- 
formly in the interval (0, 1). This value is, of course, appro- 
ximate and is in the form of an ro-digit binary fraction 
0, a (1) a (2 ) • • ■ a, m) in which each a, simulates a random 
variable with the distribution 

(11). 

v 2 2/ 

But neither this method is free of shortcomings. Firstly, 
it is difficult to check the "quality" of generated numbers. 
Checks are necessarily made periodically since the so-called 
"distribution drift" may appear due to some faults (that 
is, zeros and unities appear in some bits at nonidentical 
frequencies). Secondly, all calculations in computers are 
usually performed twice to exclude the possibility of 
occasional failure. However, it is impossible to reproduce 
the same random numbers if they are not stored in the course 
of computations. But in the case of storing them we again 
face the situation typical for tables. 

Generators of this type will undoubtedly prove iiseful 
when specialized computers will be designed for solving 
problems by means of the Monte Carlo method. But it 
would not be economical to install and operate a special 
unit in multipurpose computers in which computations, 
involving random numbers, are performed only occasional- 
ly. It is more expedient to use the so-called pseudo-random 
numbers. 

3.3. Pseudo-Random Numbers. Since the "quality" of 
random numbers used for computations is checked by means 
of special tests it is of no importance how they were genera- 
ted if only they satisfy the accepted set of tests. We may 
even calculate them by means of a prescribed formula. 
But such a formula must, of course, be a very ingenious 
one. 

The numbers calculated by means of some formula and 
simulating the values of the random variable y are called 
pseudo-random numbers. The term "simulating" means that 
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these numbers satisfy a set of tests as if they represented 
the values of this random variable. 

The first algorithm for generating random variables was 
suggested by J. Neumann. It is referred to as the mid-square 
method. Let us elucidate it by means of an example. 

Suppose we have a four-digit number y — 0.9870. Let 
us raise it to a second power. Wo get an 8-digit number 
Yo = 0.97535376. Let us select four middle digits of this 
number and assume y 1 = 0.5353. 

Now we raise y x to a second power (y] = 0.28654609) and 
again extract four middle figures. We get y„ = 0.6546. 

Further on we obtain yl = 0.42850116, yl = 0.8501; 
0.72267001, y 4 = 0.2670; y\ = 0.07128900, y 6 = 
= 0.1289, and so on.* 

However, this algorithm did not prove acceptable: the 
fraction of smaller values is higher than is necessary. Thus 
other algorithms were developed by a number of investi- 
gators. Some of them are oriented at specific features of 
particular computers. As an example we will consider one 
of such algorithms operating in the Strela computer. 

Example.** Strela is a three-address computer with a floating 
decimal point. The cell in which a number x is stored consists of 
43 binary digits (Fig. 8). The computer operates on normalized binary 



yl 



I I 1 |2|3|4| 5 



32 33 34 35|36|37 38 39 40 41 42 



v 

Mantissa 



Sign of mantissa 



Ode 



Sign of order 
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numbers x = ±<?-2 ±p where p is the order of the number and q is 
its mantissa. The ;'-th digit of the cell may contain a zero or a unity; 
we designate this value as e^. Then 



p = e 3 725 + e 38 2*- 



635 
235 



■ e 42 2°. 



* This algorithm can be written in the form y k+1 = F (y h ) where 
F designates a set of operations that must be performed over the 
number y k to obtain Yft +1 . The number y is specified. 

** See I.M. Sobol, Pseudo-random numbers for the Strela computer, 
The Probability Theory and Its Applications, 3, No. 2 (1958), 205-211. 
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The normalization condition 0.5 ^ q < 1 means that ? x is invariably 
equal to 1. The sign "+" is represented by zero, and the sign "— ", 
by unity. 

Three operations are necessary to obtain y h+1 from y h : 

1. 7/, is multiplied by a large constant, usually 10". 

2. Tho machine word representing 10 17 7 fe is shifted 7 digits to the 
left (so that the first 7 digits of the product are lost and the digits from 
36th to 42th become zeros). 

3. The absolute value of the obtained number is formed (and the 
number is normalized); this gives 7^+1- 

If we start with y = 1 , this process generates more than 80 000 
different numbers y k ; after that a cycle occurs and the numbers begin 
to repeat themselves. Various tests of the first 50 000 numbers gave 
quite satisfactory results. These numbers were repeatedly used in 
solving most diverse problems. 

Advantages of the method of pseudo-random numbers are 
fairly obvious. Firstly, only a few simple operations are 
necessary to calculate each number so that the rate of gene- 
ration is of the same order as the speed of the computer 
itself. Secondly, the program takes only a few cells of the 
storage. Thirdly, any of the numbers y k can easily be repro- 
duced. And finally, only one check of the "quality" of this 
sequence is required and thereupon this sequence can be 
used without taking any risk in calculations of similar 
problems. 

The only shortcoming of the method is that the "store" 
of pseudo-random numbers is limited. However, there are 
methods that make it possible to generate much greater 
arrays of numbers. In particular, starting numbers Yn can 
be changed. 

At present an overwhelming majority of Monte Carlo 
calculations are carried out by using pseudo-random num- 
bers. Additional information on random numbers is given 
in Sec. 10. 



Sec. 4. Transformations of Random Variables 

Solution of various problems requires simulation of various 
random variables. At an early stage of application of the 
Monte Carlo method some specialists tried to design a special 
roulette for each random variable. For instance, a roulette 
with the disk divided into unequal portions proportional 
to pi shown in Fig. 9 and functioning identically to that 
of Fig. 7, can be used to generate the values of the random 
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variable with the distribution 

/ %\ #2 %3 x k 

U.5 0.25 0.125 0.125 



(22') 



However, this proved to be absolutely unnecessary: the 
values of any random variable can be obtained by transform- 
ing a selected random variable ("standard", so to speak). 
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Usually this role is played by the random variable y, uni- 
formly distributed in (0, 1). We already know how to 
obtain values of y. 

Let us call the process of obtaining a value of some ran- 
dom variable | by means of transforming one or several 
values of y the drawing of the random variable £. 

4.1. Drawing of a Discrete Random Variable. Let us 
suppose that we want to calculate the values of the random 
variable £ with the distribution 



E = 



>.Pi 



x 2 
Pi 



Pn 



Let us consider the interval ■< y <z 1 and divide it 
into n intervals with lengths equal to p lt p 2 , • ■ ■, p n - 
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The coordinates of division points will apparently be y — 

= Pi, y = Pi + Pi, y = Pi + P2 + Pa, ■ ■ ■» y = Pi + 

+ Pa + • ■ • + Pn-i- Now we enumerate these intervals 
by the numbers 0, 1, . . ., n (Fig. 10). This completes the 
preparation procedure for drawing £. Eacli time we have 
to "run the experiment" and draw a value of £, we shall 



1 



P. + P 2 



P^'Ps 



1-Pn 



-+- 



FIG. 10 



select a value of y and fix the point y = y. If this point 
falls into the i-th interval we shall assume that £ — x t 
(in this experiment). 
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Validity of this procedure is proved very easily. Indeed, 
since the random variable y is distributed uniformly within 
(0, 1), the probability of y lying within one of the inter- 
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vals is equal to the length of this interval. Hence, 

P{0<Y<Pi} = Pi, 

P{Pi<Y<A + p 3 }--/?2, 



P {Pi + P» + • • • + Pn-l < T < !} = P<>- 

According to our procedure | = a^ if 

P1 + P2 + •••+P,-i<Y<Pi + P 2 + ■■■-!-Pn 

and the probability of this event is equal to /?;. 

Of course, Fig. 10 can be avoided when a computer is used. 
Let us suppose that the numbers x v x 2 , . . ., x n are placed 
in succession in storage cells and the probabilities p lt 
Pi + P2> Pi + Pi + P3> • • •> 1 also form a sequence in the 
data store. The block diagram of the subroutine for draw- 
ing I is shown in Fig. 11. 

Example. To draw 10 values of the random variable with the 
distribution 



( 



3 4 \ 
0.58 0.42 1 



Let us select as values of 7 ten pairs of numbers from the Table A 
on p. 74, multiplied by 0.01*. Thus, 7 -= 0.86, 0.51, 0.59, 0.07, 0.95, 
0.66, 0.15, 0.56, 0.64, 0.34. 

It is obvious that according to our scheme the value — 3 corres- 
ponds to the values of 7 smaller than 0.58 and = 4, to the values 
of y > 0.58. Hence, we obtain 0^-4,3, 4, 3, 4, 4, 3, 3, 4, 3. 

Let us note thai the order of enumerating the numbers 
x ly x 2 , . . ., x n in the distribution of | can be arbitrary, 
but it must be fixed before the drawing. 

4.2. Drawing of a Continuous Random Variable. Let us 

assume now that we need to generate the values of the ran- 
dom variable \ distributed in the interval (a, b) with the 

density p(x). 

* Two decimal digits arc sufficient in the values of y since p t 
in this example are specified only with two decimal digits. The case 
of y ~ 0.58 is possible in this approximation and should be combined 
with the case of 7 > 0.58 (because 7 = 0.00 is possible while 7 = 
= 1.00 is impossible). When multivalued 7 are used, the case of 7 = 
= p 1 is almost improbable and can be attributed to any of the ine- 
qualities. 
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Let us prove that £ values can be found from the equa- 
tion 

I 

j p(x)dx^=y, (23) 

a 

that is we must solve Eq. (23) for a selected consecutive 
value of y and calculate the consecutive value of £. 
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In order to prove it we consider the function 

X 

y= J p(x)dx. 

a 

It follows from the general properties of density, expressed 
by Eqs. (15) and (16), that 

y (a) = 0, y (b) = 1, 



and the derivative 



y' (x) =p(x)> 0. 



It means that the function y (x) monotone increases from 
to 1 (Fig. 12). Any straight line y = y where < y < 1, 
intersects the curve y = y (x) at one and only one point 
whose abscissa is taken for the value of |. Thus Eq. (23) 
always has a unique solution. 
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Now we select an arbitrary interval (a' , b') within (a, b). 
The ordinates of the curve y = y (x) satisfying the inequality 

y (a') <y<y (&') 

correspond to the points of this interval 

a' < x < b'. 




FIG. 13 

Consequently if g belongs to the interval a' < x < b', 
then 7 belongs to the interval y (a') <y <y (b'), and 
vice versa (Fig. 13). Hence 

P {a' < I < b'} = P {y (a') <y<y (&')}• 
Since 7 is uniformly distributed in (0, 1), 

6' 

V{y(a')<y<y(b')} = y(b')-y(a')= j />(*)&. 
Thus 

b' 

P{a' <£<&'}= f p(«) Ar, 



and this precisely means that the random variable £ which 
is the root of Eq. (23) has the probability density p (x). 



3* 
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Example. A random variable r\ is referred to as uniformly distribu- 
ted in the interval (a, b) if its density in this interval is constant 

p{x) = (b — a) -1 fora"0-<6. 

To draw the value of n, we apply Eq. (23): 

J b-a y - 



The integral is easily calculated: 

T) — a 



b — a 



Hence we arrive at an explicit expression for r\: 
n = a + y (b — a). 



(24) 



Other examples of application of Eq. (23) are given in 
paragraphs 5.2 and 8.3. 

4.3. Neumann's Method of Drawing a Continuous Ran- 
dom Variable. It may so happen that Eq. (23) is very hard 




FIG. 14 



to solve for |. For instance, when the integral of p(x) can- 
not be expressed via elementary functions or when the den- 
sity p(x) is specified graphically. 
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Let us assume that the random variable | is specified 
on a finite interval (a, b) and its density is limited (Fig. 14): 

p (;r) < M . 

The variable !• may be drawn as follows: 

1. Select two values y' and y" of the random variable y 
and generate a random point r(T)'; t|") with coordinates 

i\' = a + y'(b-a), v\"^y"M . 

2. If the point T lies below the curve y — p(x), assume 
| = n'; if the point T is above the curve y = p(x), reject 
the pair (y', y") and select a new pair (y' , y"). 

The validity of this method is substantiated in Sec. 9.1. 

4.4. Drawing of Normal Variables. A large number of 
various techniques exist for drawing various random varia- 
bles. We shall not consider them here. They are usually 
applied only when the methods described in Sections 4.2 
and 4.3 prove ineffective. 

Such is the case, in particular, for the normal random 
variable t,, since the equation 

— -=- \ e 2 dx — ■ y 

v -OO 

cannot be solved explicitly and the interval of possible 
values t, is infinite. 

The Table B at the end of the book (on page 74) lists 
the values (already drawn) of the normal random variable t, 
with mathematical expectation M£ = and variance 
D£ = 1. It is not hard to prove* that the random variable 

V = a + oC (25) 

will also be normal; as follows from Eqs. (10) and (11) 

MS' = a, DC = a 2 . 

Thus, formula (25) enables us to draw by using Table B 
the values of any normal variable. 

4.5. More on the Example from Sec. 1.2. Now we can 
explain how the random points for Figs. 1 and 2 were sam- 
pled. Points plotted in Fig. 1 have coordinates 

x = y', y-=y". 
* The proof is given in Sec. 9.2. 
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The values of y' and y" were calculated from groups of five 
digits of Table A: x x = 0.86515; y x = 0.90795; x 2 = 0.66155; 
y 2 = 0.66434, etc. 

It can be proved* that since the abscissae and ordinates 
of these points are independent the probability for one 
such point to fall into any region inside the square is equal 
to the area of this region. 

The points drawn in Fig. 2 have the coordinates 

*=0.5 + 0.2C, j,= 0.5 + 0.2£", 

where the values of £' and t," were selected from Table B 
in succession: 

x x = 0.5 + 0.2-0.2005, y x = 0.5 + 0.2-1.1922; 

x 2 = 0.5 + 0.2 -(-0.0077), . . . 

One of the points fell outside the square and was rejected. 
As follows from Eq. (25), the abscissae and ordinates of 
these points represent normal random variables with mean 
values a = 0.5 and variances o 2 = 0.04. 



* The proof is given in Sec. 9.3. 



CHAPTER II 

Examples of Application 

of the Monte Carlo 

Method 

Sec. 5. Simulation of a Mass Servicing System 

5.1. Description of the Problem. Consider one of the 
simplest systems of mass servicing. This system consists 
of n lines (or channels, or servicing stations) each of which 
can "serve on the customers". The system receives requests 
arriving at random moments of time. Each request arrives 
at the iVl line. If the arrival time of the /c-th request (let 
us call it T h ) finds this line free, the line starts servicing 
the request; this takes t b minutes (t b is the holding time of 
the line). If Ni line is busy at the moment T h , the request 
is immediately transferred to the N2 line. And so on ... 

Finally, if all n lines are busy at the moment T h , the 
system rejects the request. 

The problem is, what will be the (average) number of 
requests serviced by the system during the period T and 
how many rejections will be given? 

It is clear that problems of this type are encountered 
when functioning of organizations of any type, not only 
of the social service sector, is investigated. In some very 
particular cases analytical solutions were found. But in 
complicated cases (we shall discuss them below) the Monte 
Carlo method proved to be the only one available computa- 
tional method. 

5.2. Simple Flow of Requests. The first question that 
comes on when such systems are analyzed is: what is the 
flow of arriving requests? This question is answered experi- 
mentally by a sufficiently long observation of the requests. 
Investigation of the request flows in various conditions 
permitted to single out certain cases which are met with 
sufficiently frequently. 

A simple flow (or a Poisson flow) is the flow of requests 
in which the time interval t between two consecutive re- 
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quests is a random variable distributed in the interval 
(0, oo) with the density 

p(x)-=ae- ax . (26) 

The density (26) is also called the exponential distribution 
(see Fig. 15 where p(x) is plotted for a = 1 and a =' 2). 




FIG. 15 

Calculation of the expected value of i is straightforward: 

oo oo 

Mt= I xp(x)dx-~ \ xae' ax dx. 

(I 

Integration by parts (u — x, dv — ae~ ax dx) yields 



Mt- 



[ - «-«l? + j e' ax dx =[-££-] 



ax too \ 
a 



The parameter a is called the request flow density. 

The formula for drawing t is easily obtained by means 
of Eq. (23) which in our case takes the form 



1 ae ax dx- y. 
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By calculating the integral on the left-hand side we obtain 

whence 

T-, -lln(l-v). 

However, the distribution of the variable 1 — y is iden- 
tical with that of y so that instead of the above formula 
we can use the expression 

r=- -|ln Y . (27) 

5.3. The Sequence of Compulations. We shall now consi- 
der the functioning of the system of Sec. 5.1 in case of a 
simple flow of requests. 

Let each line correspond to one cell of an internal data 
store of the electronic computer; the instant when this line 
clears will be recorded in this cell. We denote the moment 
of freeing the i-th line by t t . The moment of arrival of the 
first request is chosen as the origin on the time axis: T 1 = 0. 
We can assume that all t t are equal at this moment to 7\: 
all lines are free. Computation termination time is T end = 
= T 7 ! + T. 

The first request arrives at the 7V1 line. Hence this line 
will be busy during t b . Thus, we must replace t x by a new 
value (t^ncw = T x 4 t b , add unity to the counter of ser- 
viced requests and turn to deal with the second request. 

Let us suppose that k requests have already been consider- 
ed. Now we have to draw the moment of arrival of the 
(k 4 l)-th request. We sample a consecutive value of y 
and calculate the consecutive value t = x h by means of 
Eq. (27). Then we calculate the time of arrival 

Is the first line free at this moment? To find it out we must 
check the condition 

*i<?W (28) 

If this condition is satisfied it means that at the moment 
T h+l the line is already cleared and can service the re- 
quest. We replace t x by T h+1 4- t b , add a unit to the counter 
of the serviced requests, and then switch over to the next 
request. 
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If condition (28) is not satisfied it means that at the mo- 
ment Th+i the first line is busy. We shall now check whether 
the second line is free: 

*.<ZW (29) 

If condition (29) is satisfied we replace t 2 by T k+1 + t b , 
add unity to the counter of serviced requests and switch 
over to the next request. 

If condition (29) is not satisfied, we turn to checking 
the condition 

It may happen that for all i from 1 ton 
ti > T h+i , 

i.e. all lines are busy at the moment T h + V We then add 
unity to the counter of rejections and switch over to pro- 
cessing the next request. 

Each time when Th+i is calculated we have to check the 
condition of termination of the experiment 

lh+l > Tend- 
Once this condition is satisfied, the experiment is termina- 
ted. The counter of satisfied requests and the counter of 
rejections will contain the numbers \i sat and \i re j. 

This experiment is repeated N times (each time with 
different values of y). The results of all experiments are 
averaged: 

JV 

A 



M\lre) » -jj- 2 fr 



3=1 

where \isatu) and |x rf , ; - u -> are the values of \i sa t and \i re ) 
obtained in the j'-th experiment. 

Fig. 16 shows the block-diagram of a program realizing 
these computations (if necessary, we can retrieve the values 
l*>sat<j> and \x r ej(j) of individual trials in the "end of expe- 
riment" block). 

5.4. More Complicated Problems. It is not hard to show 
that the same method permits computation for much more 
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complicated systems. For instance, the holding time t b 
may be not constant but random and vary from line to line 
(this corresponds to different types of equipment and to 
unequal skills of the staff). The scheme of computations 
remains basically unchanged but the values of t h will have 
to be drawn in each trial and the formula for drawing will 
be specific for each line. 

We can also consider the so-called systems with waiting 
in which the rejection is not immediate: -the request is 
kept for some time t r (time the request remains in the system) 
and is serviced if one of the lines is cleared during this 
time. 

Another type of systems, in which the next request is 
serviced by the first line to be cleared, can be considered. 
We can take into account a random failure of each individual 
line and a random time interval necessary to repair it. 
We can take into account time variations of the request 
flow density, as well as many additional factors. 

As a matter of course, a certain price has to be paid for 
simulation of such systems. Practically valuable results 
are only obtained if a good model is selected. For this pur- 
pose we have to study thoroughly actual request flows, to 
carry out timing of functioning of individual equipment 
units, etc. 

Generally speaking, we have to know the probabilistic 
laws governing the functioning of separate parts of the 
system. Then the Monte Carlo method makes it possible to 
calculate the probabilistic laws of the entire system ope- 
ration, regardless of its complexity. 

Such methods of calculation are extremely helpful in 
designing economic units: expensive (and sometimes 
even impossible) actual experiment is substituted by a 
simulated computer experiment in which various versions 
of work organization or equipment operation are simulated. 

Sec. 6. Computation of Quality 
and Reliability of Complex Devices 

6.1. The Simplest Scheme of Computing Quality. Let us 

consider a device S consisting of a certain (possibly large) 
number of elements. If, for instance, S is an electric appara- 
tus, then its elements can be resistors (B (k) ), capacitors 
(C<fc>)> electron tubes, etc. Let us assume that the quality 
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of the device is determined by the value of a single output 
parameter V which can be calculated from the known para- 
meters of all the elements 



V :i -~ f (-ft(D) -ftcSl! • • ■> ^<1>> ^<2>i 



•)■ 



(30) 



If, for example, V is the voltage across an operating part 
of an electric circuit, the equations for this circuit can be 
written by using Ohm's laws so that V is obtained by solv- 
ing these equations. 




FIG. 17 



In practice, however, the parameters of the elements are 
not exactly equal to the indicated values. For instance, 
the resistance of the resistor shown in Fig. 17 may vary 
within the range from 20.9 to 23.1 kQ. 

The question now is: how will V be affected by deviations 
of the parameters of all elements from their rated values? 

We can try and estimate the limits of V fluctuations by 
choosing the "worst" values of the parameters for all ele- 
ments. However, it is not always known which set of para- 
meters is the "worst"? Moreover, if the total number of ele- 
ments is large, such an estimate will be much too high: 
it is indeed very unlikely that all parameters will take on 
their worst values simultaneously. 

It is, therefore, more reasonable to consider the parame- 
ters of all elements and the variable V itself random and 
to try and estimate the expected value MV and the vari- 
ance DV. The value MV is the mean value of V for the whole 
lot of items and DV shows what deviations of V from MV 
will be encountered in practice. 

Let us recall (it was mentioned in Sec. 2.2) that 



MF^/(Mi? (1) , M/? (2) , ...; MC (1) , MC (2) , 



•)• 



The distribution of V cannot be computed analytically 
if the function / is complicated to even the smallest extent. 
Sometimes we can do it experimentally by scanning a large 
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batch of final items of the device. However, this is not always 
possible; and it is never possible at the design stage. 

Now let us attempt to apply the Monte Carlo method.. 
For this we need: (a) knowledge of the probabilistic charac- 
teristics of all elements; (b) knowledge of the function / 
(rather, we must be able to calculate V values from all 
specified values R (1) , i?< 2 >, • • •', C a> , C (2) . . .; . . .). 

Probabilistic distributions of parameters for each indivi- 
dual element can be obtained experimentally, by scanning 
a large batch of such elements. These distributions are very 
often normal. For this reason many investigators use the 
following approach: they consider, for example, the resi- 
stance of the element shown in Fig. 17, as a normal random 
variable p with the expected value Mp = 22 and 3a = 1.1 
(it will be recalled that it is practically improbable that 
we obtain in one trial the value of p that deviates from Mp 
by more than 3a; see Eq. (20)). 

The scheme for the computation is quite simple: first 
the value of the parameter is drawn for each element; then 
the V value is calculated by formula (30). Having repeated 
this experiment N times we obtain values V lf V 2 , ■ ■ ., V N 
and can approximately assume that 

N 

N Zl " 
i=i 

j=l 3=1 

At large values of N the factor \I(N — 1) in the last 
formula can be replaced by 1A/V; then this formula will di- 
rectly follow from Eqs. (8) and (9). Mathematical statistics 
proves that the factor l/(N — 1) provides better results 
for small N. 

6.2. Examples of Calculating Reliability. Suppose we 
want to estimate the mean duration of failure-proof opera- 
tion of the apparatus provided all characteristics of the 
failure-proof operation of each of the elements are known. 

If we assume that the duration of failure-proof operation of 
each element t ik) is fixed, then the calculation of time t of 
failure-proof operation of the whole unit presents no pro- 
blems. For example, failure of any one of the elements of 
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the unit drawn in Fig. 18 entails the failure of the whole 
unit, i.e. 

t^--min{t (V ; i (2) ; t i3) ; i (4) ). (31) 

For the unit shown in Fig. 19 where one of the elements has 
a stand-by element 



t= min [t (i) , t i2) ; max(£ (3) ; t a) ; i (5; ], 



(32) 



since when one element, iV3 for example, will break down, 
the unit will continue functioning on the remaining ele- 
ment iV4. 
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FIG. 18 



In reality the duration of failure-proof operation of any 
single element is a random variable (Jl) . When we say that 
the service life of an electric bulb is 1000 hours it is only 




FIG. 19 



the mean value M0 of the variable 0: everyone is aware 
that one bulb burns out sooner while another one (exactly 
identical with the former) lives longer. 

When distribution densities of Q ik) are known for each 
of the elements of the unit, M9 can be computed exactly 
as it was performed in Sec. 6.1. Indeed, (h) can be drawn 
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for each element; let it be t (h) . Then the value t of the ran- 
dom variable can be calculated by means of relevant for- 
mulas (31) or (32). Having rerun this trial a sufficient number 
of times (TV) we can accept that 



M6«^2^ 



j= 



where tj is the value of t obtained in the ;-th trial. 

It should be mentioned that the problem of distribution 
of service life 0< ft) of individual elements is not so simple 
as one might think: running of an experiment involving 
maximum longevity elements is cumbersome because we 
need to wait until a sufficiently large number of elements 
break down. 

6.3. Further Possibilities of the Method. The examples 
given above demonstrate that in principle the techniques 
of computing quality of the articles being designed is quite 
simple. One must know probabilistic characteristics of all 
elements of the unit and be able to calculate the value of 
interest as a function of the parameters of these elements. 
Then the random character of these parameters is taken 
into account by simulation. 

Simulation can yield much more useful information than 
just the expected value and the variance of the variable 
of interest. Suppose, for example, that we have sampled a 
large number N of values U x , U 2 , . . ., U N of a random 
variable U. On the basis of these numbers we can plot an 
approximate distribution density of U. In fact this is the 
domain of statistics because we touched upon processing 
of experimental data (produced in computer-simulated 
experiments). Thus we will restrict ourselves to the presen- 
tation of one specific example. 

Let us assume that all in all we obtained N = 120 values 
U u U 2 , • • ., U 120 of the random variable U, all distribu- 
ted in the range 

1 < U } < 6.5. 

Let us divide the interval 1 < x <C 6.5 into 11 (or any 
other number, not too large and not too small) equal inter- 
vals with lengths Ax = 0.5 and count how many values 
of Uj fall into each of the intervals. Corresponding numbers 
are indicated in Fig. 20. 
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The frequency of falling into any one of the intervals 
is calculated by dividing these numbers by N = 120. In 
our case the frequencies are: 0.017; 0; 0.008; 0.12; 0.20; 
0.27; 0.14; 0.16; 0.06; 0.008; 0.017. 
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FIG. 20 



Now we draw a rectangle over each of the intervals with 
the area equal to the frequency corresponding to this inter- 
val (Fig. 21). In other words, the height of each rectangle 
is equal to frequency divided by Ax. The step line obtained 
is called a histogram. 




The histogram is an approximation to the unknown den- 
sity of the random variable U. Thus, for instance, the area 
of this histogram between x = 2.5 and x = 5.5 yields an 



4-01015 
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approximate value of the probability 

P{2.5 < £/<5.5} « 0.94. 

Hence, we can assume on the basis of the above calculation 
tliat with the probability approximately equal to 0.94 
the value of U lies within the interval 2.5 < £/<<5.5. 

For the purpose of comparison the density of a normal 
random variable £' with parameters a = 3.85, o = 0.88 
is plotted in Fig. 21. If the probability of £' being inside 
the interval 2.5 < £' < 5.5 is calculated from the normal 
density, we obtain a fairly close value 0.91*. 

6.4. Note. Unfortunately, calculations of this type are 
still fairly scarce. It is not easy to say what is the main 
reason for this. The most likely explanation is that design- 
ers simply are not aware of this possibility. 

Besides, before performing computations concerning de- 
vices of any type we have to analyze the probabilistic cha- 
racteristics of all elements incorporated into these devices. 
This means a lot of work. But once these characteristics 
are known, we can evaluate the quality of any device con- 
sisting of these elements. The variation of the quality caused 
by replacing some of the elements by other ones can be 
estimated as well. 

* The course of calculations. According to Eq. (14) we have 

5.5 (.v-q)2 

P{2.5<£'<5.5} = i —=- [ e 2o2 dx. 

ay 2n 2 J 5 

We change the variable in the integral x — a = at. This yields 

— ■ l f 



P{2.5<£'<5.5} = — ~ [ e z dt, 

y 2n j 



n 

where t^ — — ■•■-■ —1.54, t 2 = — — 1.88. The last integral 

is evaluated by using one of the tables of the so-called error function 
which lists the values of the function 



X _J2_ 
Cl)(i)- — ^=- \ e - dt. 





x 

J 2n J 



We obtain 

P {2.5 < Z.' < 5.5} --.= 0.5 [<I> (1.88) + 0> (1.54)] = 0.91. 
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We should hope that in the nearest future such calcula- 
tions will become a routine business, and the probabilistic 
characteristics of the elements will invariably be supplied 
by their manufacturers. 

Sec. 7. Computation of Neutron Transmission 
Through a Plate 

The probabilistic laws of interaction of an individual 
elementary particle (neutron, photon, meson, etc.) with 
matter are known. Usually it is needed to find the macro- 
scopic characteristics of processes in which the number of 
participating particles is enormous, such as densities, fluxes 
and so on. This situation is quite similar to that encountered 
in Sees. 5 and 6, and offers a very convenient object for the 
Monte Carlo calculations. 

The neutron physics is probably the field where the Monte 
Carlo method is used more frequently than anywhere else. 
We shall consider the simplest version of the problem of 
neutron transmission through a plate. 

7.1. Statement of the Problem. Let the flux of neutrons 
with energy E be incident on a homogeneous infinite plate 
^ x ^ h. The angle of incidence is 90°. Upon collision 
with atoms of the plate material neutrons may be either 
elastically scattered or captured. Let us assume for simpli- 
city that energy of neutrons is not changed in scattering 
and that any direction of "recoil" of a neutron from an atom 
is equally probable (this is sometimes the case in neutron 
collisions with heavy atoms). Fig. 22 illustrates possible 
fates of a neutron: neutron (a) passes through the plate, 
neutron (b) is captured, and neutron (c) is reflected by the 
plate. 

It is necessary to calculate the probability of neutron 
transmission through the plate p + , the probability of neu- 
tron reflection by the plate p- and the probability of neu- 
tron capture inside the plate p°. 

Interaction of neutrons with matter is characterized in 
the case under consideration by two constants 2c and 2s 
which are called the capture cross-section and the scattering 
cross-section, respectively. 

The sum of these two cross-sections is called the total 
cross-section 
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The physical meaning of these cross-sections is as follows: 
the probability of neutron capture upon collision with an 
atom is equal to ^</2j, and the probability of scattering 
is equal to 2j S /2j- 




FIG. 22 

The free path length of a neutron A (i.e. the length of the 
path from one collision to another) is a random variable. 
It can take on any positive value with the probability 
density 

It is easily seen that this density of A. coincides with the 
density (26) of the random variable t of a simple request 
flow. Similarly to Sec. 5.2 we can immediately put down 
an expression for the mean free path length 

MA, = 1/2 

and the formula for drawing A, values 

X= -4r lny - 
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We now have to find out how to select a random direction 
of a neutron after scattering. As the problem is symmetrical 
about the x-axis, the direction is completely denned by an 
angle cp between the neutron velocity and the Ox-axis. 
It can be proved* that the requirement of equal probability 
of all directions is equivalent in this case to the require- 
ment that the cosine of this angle u. = cos cp be uniformly 
distributed over the interval ( — 1, 1). Formula (24) for 
a — — 1, 6 = 1 yields the expression for drawing \i values: 

u. = 2 Y - 1. 

7.2. Computation by Means of Simulation of Real Trajec- 
tories. Let us assume that a neutron has undergone the i-th 




FIG. 23 

scattering inside the plate at the point with the abscissa 
x h and started moving after it in the direction \i lt . 
We draw the free path length 

X h = -(1/2) In y 

and calculate the abscissa of the next collision (Fig. 23) 

x hH '- x k~\~ ^k n h- 

* The proof is given in Sec. 9.4. 
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We check the condition of neutron transmission through 
the plate: 

If this condition is satisfied, the computation of the trajecto- 
ry is terminated and unity is added to the counter of trans- 
mitted particles. Otherwise \vc check the condition of reflec- 
tion: 

x h+i < 0. 

If this condition is satisfied, the computation of the trajecto- 
ry is also terminated and unity is added to the counter of 
reflected particles. If this condition is not satisfied, i.e. 
^ x h+1 ^ h, it means that the neutron has undergone 
the (k -j- l)-th collision inside the plate and we have to draw 
its "fate" in this collision. 

According to Sec. 4.1 we choose the next value of y and 
check the capture condition 

y < Se/2- 

If this last inequality is satisfied, the trajectory is terminated 
and we add unity to the counter of captured particles. 
Otherwise we assume that the neutron had been scattered 
at the point with the abscissa x h+1 . In this case we draw a new 
direction of the neutron velocity 

Hfe+i = 2Y-l 

and then repeat the whole cycle (but, of course, with new 
values of y). 

In all of the above formulas we added no subscript to y 
because it is meant that each one value of y is used only once. 
Three values of y are necessary to compute one lap of the 
trajectory. 

The initial values for each trajectory are: 

*o =0, Ho = 1. 

After N trajectories are sampled we obtain that N + 
neutrons were transmitted through the plate, N~ neutrons 
were reflected and N° neutrons captured by it. It is obvious 
that the probabilities we seek are appoximately equal to the 
ratios 
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Fig. 24 shows the block diagram of the program of compu- 
tations for this problem. The subscript / is the trajectory 
number and the subscript k is the collision number (along 
the trajectory). 

This computation procedure, though quite natural, is 
far from being perfect. In particular, by applying this 
method we run into difficulties if the probability p + is very 
small. However, this is just the case that we have in design- 
ing a radiation-proof shielding. 

There exist some more "ingenious" versions of the Monte 
Carlo method which make it possible to carry out calcula- 
tions in this case as well. We shall now briefly discuss one 
of the simplest versions of computations involving the so- 
called "weights". 

7.3. Computation Scheme Using Weights which Replace 
Absorption of Neutrons. Let us consider now the same pro- 
blem of neutron transmission. Let us suppose that a "packet" 
consisting of a large number w of identical neutrons moves 
along a single trajectory. The average number of neutrons 
captured from the "packet" in a collision at the point with 
the abscissa x { , will be equal to w (E c /2), and the average 
number of scattered neutrons, to w (S,/J.). 

We add the quantity w (2 c /2) to the counter of cap- 
tured particles and will follow the scattered "packet" under the 
assumption that the remaining "packet" was scattered as 
a whole in one direction. 

All computation formulas given in Sec. 7.2 remain un- 
changed but the number of neutrons in the "packet" upon 
each collision will be reduced to 

Wft + l = H.'ft(2s/S)» 

since part of the "packet", containing w k (2 c /2) neutrons 
will be captured. The trajectory cannot be terminated by 
capture any more. 

The quantity w h is usually referred to as the weight of 
a neutron; and rather than speak about a "packet" consisting 
of w h neutrons we speak about one neutron with the weight 
w h . The initial weight w is usually assumed to be 1. This 
is not contrary to what was said about a "large packet" be- 
cause it is easily seen that all w k obtained in computations 
of a single trajectory contain w^ as a common multiplier. 

The block diagram of a program realizing this computa- 
tion is shown in Fig. 25. This diagram is not a shade more 
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complicated than thai shown in Fig. 24. However, it can be 
proved* that the computation of p + by the method of this 
section is always more expedient than by the method covered 
in Sec. 7.2. 

7.4. Note. The number of various methods of calculations 
using different weights is quite considerable. We cannot 
dwell on all of them here. 

We will only emphasize that the Monte Carlo method 
makes it possible to solve much more complicated problems 
concerning elementary particles: the investigated medium 
may consist of different substances and be of an arbitrary 
geometric shape; the energy of a particle may change in 
each collision. Many other nuclear processes may be taken 
into account (for example, the possibility of atom fission 
upon collision with a neutron and subsequent formation of 
additional neutrons). Conditions at which chain reaction 
is initiated and maintained can also be calculated, etc. 

Sec. 8. Calculation of a Definite Integral 

The problems considered in Sees. 5, 6 and 7 were of pro- 
babilistic nature so that the application of Monte Carlo 
techniques to their solution appeared fairly natural. Here we 
consider a standard mathematical problem: approximate 
computation of a definite integral. 

Since the calculation of definite integrals is equivalent 
to the calculation of areas, we might utilize the method 
presented in Sec. 1.2. As alternative to this method we 
shall describe another more efficient method making possible 
the construction of various probabilistic models to solve 
this problem by the Monte Carlo method, and indicate how 
to select a "better" model among these. 

8.1. Method of Calculations. Let us consider a function 
g (x) defined on the interval a < x << b. We have to appro- 
ximate the integral 

l=^g{x)dx. (33) 

a 

Let us choose an arbitrary distribution density pi (x) 
specified on the interval (a, b) (i.e. an arbitrary function 
Pl (x) satisfying the conditions (15) and (16)). 

* The proof is given in Sec. 9.5. 
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Along with the random variable | defined on the interval 
(a, b) with the density Pi (x), we shall need another random 
variable 

i = * (IVpid). 

According to (18) 

b 

MT| = j [g (x)/pt (*)] P-, (x) dx ---- 1. 

a 

Let us consider now N identical independent random varia- 
bles Th, ti 2 , . . ., r\ N and apply the central limit theorem 
of Sec. 2.4 to their sum. In this case Eq. (21) will take the 
form 

''{|4-2 t u- / |< 3 /¥}~ - 997 - ( 34 ) 

The last relation means that if we sample N values | lt 
E2. • • •» lv> then for sufficiently large N 

T^«'- (35) 

It also shows that the error of approximation (35) will not 
exceed, with very high probability, the value 



3 /DnA/V. 

8.2. How to Choose a Computation Scheme. Wo have 
seen that any random variable |, defined in the interval 
(a, b) can be used to calculate integral (33). In all cases 

Mt) = M [g (l)/ Pl (I)] = /. 

But the variance Dn, and with it the estimated error of 
(35) depend on what specific variable £ is used. Indeed 

b 

Dn -- M (ti 2 ) - P - ( [g 2 (x)/p l (x)\ dx- P. 

a 

It can be proved* that this expression takes on the mini- 
mum when pi (x) is proportional to | g (x) | . 

Of course, too complex pg (x) should not be chosen because 
the procedure of drawing | values will become excessively 

* The proof is given in Sec. 9.6. 
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time-consuming. But the suggestion given above is a useful 
guide-line in the selection of p\(x) (see the Example of 
Sec. 8.3). 

In practice integrals of type (33) are not calculated by 
means of the Monte Carlo method: more accurate methods, 
namely quadrature formulas are available. But the situa- 
tion changes drastically when we consider multiple inte- 
grals: quadrature formulas become extremely complicated 
while the Monte Carlo method remains almost unchanged. 

8.3. Numerical Example. Let us approximately calculate the 
integral 

It/2 



/= \ sin x dx. 

The exact value of this integral is known 

it/2 



\ sin x ii = [ — cos;r]* /2 = l. 



We shall use two different random variables % for this calculation: 
one with constant density p^(x) = 2/n (i.e. £ is uniformly distri- 

8x 




FIG. 2D 



buted in the interval (0, jt/2) and the other with linear density 
Pg (x) — 8j/ji 2 . Both those densities and the integrand sin x are plotted 
in Fig. 26. This figure shows that the linear density agrees better with 
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the suggestion of Sec. 8.2 that proportionality of p% (x) and sin x be 
desirable. Thus we can expect that the second way of computation will 
produce a better result. 

(a) Let p| (x) = 2/jt in the interval (0, n/2). The formula for 
drawing | can be obtained from Eq. (24) with a — and b — n/2: 

1 = jiy/2. 
Now formula (35) takes the form 

N 
7=1 

Let us set N = 10. For the values of y we use triads of numbers 
from Table A (multiplied by 0.001). Intermediate results are listed 
in Table 1. 

TABLE i 



j 


1 


2 


3 


4 


5 


Vj 


0.865 
1.359 

0.978 


0.159 
0.250 
0.247 


0.079 
0.124 
0.124 


0.566 
0.889 
0.776 


0.155 
0.243 
0.241 


/ 


6 


7 


8 


9 


10 


yj 


0.664 
1.043 
0.864 


0.345 
0.542 
0.516 


0.655 
1.029 
0.857 


0.812 
1.275 
0.957 


0.332 
0.521 
0.498 



The final result is 

/ x 0.952. 

(b) Now we set p§ (x) = 8z/.n 2 . To draw g we use Eq. (23) 
| 

\ (8x/n i )dx = y, 
o 
whence, after simple calculations, we obtain 

Eq. (35) takes now the following form 

N 



t JL_ V sin S; 

" 8AT Zj E, 



;=1 
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Let N — 10. We choose the same values of y as in the case (a). 
The intermediate results are listed in Table 2. 



TABLE 



i 


1 


2 


3 


4 


5 


v; 


0.865 


0.159 


0.079 


0.566 


0.155 


\i 


1.461 


0.626 


0.442 


1.182 


0.618 


sin %j 


0.680 


0.936 


0.968 


0.783 


0.937 


J 


6 


7 


8 


9 


10 


yj 


0.664 


0.345 


0.655 


0.812 


0.332 


li 


1.280 


0.923 


1.271 


1.415 


0.905 


sin| 7 - 


0.748 


0.863 


0.751 


0.698 


0.868 



The final result is 

J x 1.016. 

As we have expected the second way of computation yielded a 
more accurate result. 

8.4. On the Error Evaluation. As noted in Sec. 8.2 the 
absolute error of computation of the integral / by means of 
formula (35) cannot practically exceed the value 3 Y Dr\/N. 
But the actual error is, as a rule, appreciably smaller than 
this estimate. Thus in practice another quantity is used 
to characterize the error, namely, the so-called probable 
error 

o pro6 - 0.675 YDrjjN. 

The actual absolute error depends on random numbers 
used in the computation and may prove to be 2-3 times 
greater or several times smaller than 8 pTO b- Thus, 6 Pro( , 
estimates not the upper limit of the error but only its order 
of magnitude*. 

Let us go back to the example of Sec. 8.3. The values listed in 
Tables 1 and 2 can be used for an approximate evaluation of variances 



* Definition of the probable error is given in Sec. 9.7. 

62 



Dt| {or both methods of computation.] The relevant computation 
formula is given in Sec. 6.1.* 

The approximate variances Dn, probable errors 8 Prob calculated 
from these variances, and actual errors of calculations 6 cn / c are listed 
in Table 3. 



TABLE 3 



Method 


Dri 


6 prob 


6 calc 


(a) 
(b) 


0.256 
0.016 


0.103 
0.027 


0.048 
0.016 



It can be seen that 6 CQ ; C are indeed of the same order as 6 prob . 



* In the case (a): 

10 10 

D n-^[S( sln ^ 2 -4r(S sin ^) 2 j = 



= %- (4.604 - 3.670) = 0.256. 



In the case (b) 



sing,- \2- 



n n* r V / s '"5j \ 2 1 (y-. sinjj ^2-, 

D ^9^4-|.2i(-|-) — io(S— IT) J 



=---—- (6.875 — 6.777) =.- 0.016. 
57b 
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Appendix 



Sec. 9. Proofs of Selected Statements 

This section contains proofs of certain statements of- 
fered in the foregoing sections, these proofs either appear 
somewhat bulky for a popular text or expect a more exten- 
sive knowledge of the probability theory. 

9.1. Proof of von Neumann's Method of Drawing a Ran- 
dom Variables (Sec. 4.3). The random point T is uniformly 
distributed in the rectangle abed (Fig. 27) whose area is 
M (b — a)*. The probability for the point T to fall below 
the curve y = p (x) and thus not to be rejected is equal to the 
ratio of the areas 

b 

\ p (x) dx 



M (b — a) ' M (b—a) ' 



The probability for the point T to fall below the curve 
y = p (x) in the interval a' <C x < b' is also equal to the 
ratio of areas 

b' 

j p (x) dx 



M (b-a) ' 

Hence, the fraction of g values comprised within the 
interval (a', b') among all selected values of \, is equal 
to the ratio 

6' 

J p[x) dx b' 

M (b-a) : M (b-a) = J P ^ dx ' 
a' 

that is what is required to prove. 
* Cf. Sec. 9.3. 
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9.2. Distribution Density of the Variable g' = a + ag 
(Sec. 4.4). It is assumed that the variable £ is normal and 




has the expected value M£ = and the variance D£ = d, 
so that its density is 

To calculate the distribution density of the variable £' 
we select two arbitrary numbers x 1 < z 2 and calculate the 
probability 

P {*i < C < «2> = P {a?i < a + o£ < z 2 } = 

_p / x t — a v :k 2 — a •> 



Hence 



5-01015 



(.r 2 -a)/o ^^2_ 

P{x,<g'<af 2 } = -7==- f e" 2 <fe. 

1/ 2n J 



(x,-o)/a 



05 



Let us change variable x for x' = a -f- ox in the last inte- 
gral. We obtain 

Y{x,<l!<x^^—^r-\e ^ dx', 

x i 

whence stems (cf. Eq. (14)) the normality of the variable £' 
with the parameters M£' = a, D£' = o 2 . 

9.3. Uniform Distribution of Points (y , Y ") Inside a 
Square (Sec. 4.5). Since the coordinates of the point (y', y") 
are independent, its density p (x, y) is equal to the product 
of densities 

p(x, y) = p Y (x)p r (y). 

Each of these densities is identically equal to 1. Henee, 
p (x, y) = 1 (for < x < 1 and < y < 1). This is exactly 
what we define as uniform distribution of the point (y' , y") 
in the unit square. 

9.4. Selection of a Random Direction (Sec. 7.1). Let us 
specify the direction by a unit vector emerging from the 
origin of coordinates. The ends of these vectors lie on the 
surface of the unit sphere. The expression "any direction 
is equally probable" means that the end of the vector is a 
random point Q uniformly distributed over the surface of 
this sphere: the probability that Q will be found in any ele- 

JO 

mont of the surface dS is equal to -r- • 

Let us choose spherical coordinates (cp, ty) with the 
polar axis Ox on the surface of the sphere (Fig. 28). Then 

dS = sin cp-dcp-dij), (36) 

where ^C 9 ^ n, 0^ip< 2jt. 

Let us denote by p (<p, op) the density of the random point 
(cp, \p). It follows from Eq. (36) and from the requirement 

P(% y)d(p-d\p = ^ 

that 

P(«P. 1-)=-^- (37) 

(G 



Densities of (p and \p can be easily calculated from the joint 
density of these variables: 

2jt 

M<P)= jp( f P, *)d* = i x L : (38) 







p*!(*)=j p(<p> <p)<*<p = -Jp 



(39) 



The equality p (cp, ip) = ^ (cp) p^, (ap) demonstrates that (p 
and o|) are independent. 




FIGX28 

It is obvious that ip is uniformly distributed in the inter- 
val (0, 2n) and the formula for drawing ip should be writ- 
ten in the form) 

i|j = 2ny. (40) 

The formula for drawing (p will be found by means of 
Eq. (23): 

<p 

y I sin x dx — y, 



5* 
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whence 

cos (p = 1 — 2v. (41) 

Eqs. (40) and (41) enable us to sample (draw) a random 
direction. Of course, the values of y in these formulas must 
be different. 

Equation (41) differs from the last formula in Sec. 7.1 
only in that 1 — y is replaced for y, both variables having 
identical distributions. 

9.5. Advantages of the Computation Method Using 
Weights (Sec. 7.3). Let us introduce random variables v 
and v' equal to the number (weight) of neutrons transmitted 
through the plate and obtained from one trajectory computed 
by the methods of Sec. 7.2 and 7.3, respectively. 

It is physically obvious that 

Mv = Mv' = p + . 

(This statement is rigorously proved in the author's mono- 
graph [3].) 

Since v can take on only two values, or 1, the distribu- 
tion of v is specified by the table 

Taking into account that v 2 = v, we easily find that 
Dv = p + — (p + ) z . 

It is clear that the variable v' can take on an infinite set 
of values: w = 1, w i , w 2 , . ■ ., w h , . . ., as well as 0. Thus 
its distribution is specified by the table 



v = 



The values of q ( are of no interest for us because the variance 
can be written in all cases as 



w 


Wi 


w 2 , 


. . w h . 


.. 


?0 


Qi 


qi ■ 


. . q k .. 


• • q. 



Dv'= 2 wlq h -tf)\ 

oo 

We note that all w h < 1 and that 2 w hqu = Mv' = p + , 

and thus obtain the inequality Dv '^ p* — (p + ) 2 = Dv. 
The fact that the variance of v' is always smaller than 
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that of v demonstrates that the method of Sec. 7.3 is always 
better for computation of p + than that of Sec. 7.2. 

The same conclusion is valid as regards the calculation 
of p~, and, provided the absorption is not too high, as 
regards the calculation of p°. 

9.6. The Choice of the Best | (Sec. 8.2). An expression 
for the dispersion Dn was obtained in Sec. 8.2. To find 
the minimum of this expression among all possible choices 
of pi (x) we will use a well known inequality 

b b b 

| 1 | u (x) v (x) | dx\ <; \ u 2 (x) dx- \ v 2 (x) dx. 

a' a a 

Let us set u = g (x)/Y p% (x) and v = V Pi ( x )> th en tms 
inequality yields 

[l\g{')\dxy<]^dx]p,(x)d,=\^dx. 

o, a a a 

Hence 

b 

Dti>[j | g(s) | <**]*-/*. (42) 

a 

It now remains to prove that the lower boundary is taken 
on when p| (x) is proportional to | g (x) \ . 
Let 

Pt(*)= >" WI • (43) 

\ \g(x)\ dx 
a 

It can be easily calculated that for this density p$ (x) 

i[^]*-[ii'wi""]'- 

a a 

and the variance Dn is indeed equal to the right-hand side 
of the expression (42). 

Note that it is practically impossible to select the "best" 
density (43) for these calculations: this would require that 

6 

we know the value of the integral \ | g (x) | dx. However, the 

a 

evaluation of the last integral is a problem equivalent to 

09 





that of evaluation of I g (x) dx. Thus we have restricted 

a 

ourselves to the suggestion given in Sec. 8.2. 

9.7. Definition of the Probable Error (Sec. 8.4). Let £ 
be the normal random variable defined in Sec. 2.3. It is 
easily calculated that whatever a and o we have for r = 
= 0.675a 

a + r 

\ Pz(x) dx ----- 0.5. 

a-r 

Hence 

P{|E-a|<r} = P{|C-a|>r} = 0.5, 

i.e. deviations exceeding r and those smaller than r are 
equally probable. Thus, the characteristic r is referred 
to as the probable error of the random variable £. 

In Sec. 8.1 we calculate an approximately normal variable 
p = (1/JV) (r\ 1 + r\ 2 -f . . . + 11 jv)- Its expected value is 
a = Mp = /, and the variance is a 2 ss Dp = Dn/iV. Hence, 
the probable error of the variable p is approximately equal 
to 0.675 ^DV7V7 

Sec. 10. On Pseudo-Random Numbers 

Most of the algorithms for generation of pseudo-random 
numbers have the form 

y k+i -F(y h ). (44) 

If the starting number y is fixed, all the successive numbers 
Tii ?2» • • • are calculated by the same formula (44) for 
k — 0, 1, 2, ... . Both algorithms considered in Sec. 3.3 
also have form (44) although a set of operations, which are 
to be carried out over the argument x to obtain a value of y, 
was given instead of specifying the function y = F (x) 
in an analytical form. 

10.1. What Should Be the Function F (x). The example 
below helps to elucidate the nature of one of the basic dif- 
ficulties of selecting F (x). 

Example. Let us prove that the function F (x) plotted in Fig. 29 
cannot be used to generate pseudo-random numbers by means of 
formula (44). 

70 



Indeed, let ns consider the points with Cartesian coordinates 
(Vi. 7a). (7a. 74). (75- Ye). • • • 

inside the unit square < x < 1, < 1/ < 1. Since here we have 
Y2 ~ F (7i). 7i = F (73). 7b = F (y s ), . . ., all these points arc located 
on the curve y = F (,r). This is of course unacceptable because true 
random points 'must uniformly fill the whole square (see Sec. 9.3). 




fig. 29 

The above example leads to the conclusion that the func- 
tion y = F (x) can be successfully used in formula (44) 
only if its graph provides sufficiently dense filling of the 
whole square! 

An example of a function possessing this properly is 

y = {^}, (45) 

where g is a very large number and {z} is the fractional part 
of the number z. Function (45) is plotted in Fig. 30 for g — 
= 20. The reader can imagine what does this plot look like 
at g = 5 17 . 

10.2. The Congruential Method (Method of Residues). 
The most widespread algorithm for generation of pseudo- 
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random numbers was suggested by D. Lehmer. This algo- 
rithm is based on function (45) but for the sake of conve- 
nient realization in computers the calculation formulae are 
written somewhat differently. 



v,\ 



'={20x} 



FIG. 30 

A sequence of integer numbers m h is generated in which 
the starting integer m = 1 is fixed and all subsequent inte- 
gers m u m 2 , . . . are calculated by means of one formula 

m h+1 = 5"m k (mod 2 40 ), (46) 

for k = 0, 1, 2, . . .; from the numbers m/, we calculate 
pseudo-random numbers 

Y/.-2-*°m fc . (47) 

Formula (46) means that the number m J+ i is equal to the 
residue obtained as a result of dividing 5 17 /n ft by 2 40 . In the 
theory of congruence (see any textbook on the theory of 
numbers) this residue is referred to as the least positive 
residue modulo 2 40 , Thus both terms for this algorithm— 
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congruential method and method of residues — originate from 
this theory. 

Both formulas (46) and (47) are very easily realized in 
electronic computers, operating with 40-bit numbers, 
by means of the multiplication command with double 
precision: low-order digits of the product must be used*. 

The period of the sequence m k is equal to 2 38 . 

In computers operating with 36-bit numbers we use 5 15 
and 2 36 instead of 5 17 and 2 40 in formulas (46) and (47); 
the period of the sequence m h is in this case 2 34 . 



* The computer Strela (see Sec. 3.3) has no such a command. 



RANDOM NUMBER TABLE 



Table A. 400 random numbers * 



86 515 


90 795 


66 155 


66 434 


69 186 


03 393 


42 502 


99 224 


41686 


42163 


85181 


38 967 


86 522 


47171 


88 059 


89 342 


72 587 


93 000 


89 688 


78416 


52 452 


42 499 


33 346 


83 935 


76 773 


97 520 


27 256 


66 447 


04825 


82134 


80 317 


75 120 


87 113 


84 778 


45 863 


24 520 


84 754 


57 616 


38132 


64 294 



Table B. 88 random normal variables ** 



0.2005 


1.1922 


—0.0077 


0.0348 


1.1609 


—0.6690 


-1.5893 


0.5810 


0.5864 


—0.9245 


0.0904 


1.5068 


0. 1425 


—0.2863 


1.2809 


0.4043 


0.9516 


—1.7708 


2.8854 


0.4686 


0.5863 


0.8574 


-0.5557 


0.8115 


1.1572 


0.9990 


—0. 1032 


0.5405 


1.4428 


—0.5564 


—0.5098 


—1.1929 


0.3924 


1.7981 


0.6141 


—1.3596 


0.8319 


0.4270 


—0.8888 


0.4167 


0.9780 


-0.7679 


0. 8960 


0.5154 



* Random numbers simulate the values of a random variable with 
** Normal variables simulate the values of a normal (Gaussian) random 



56 558 


12 332 


94 377 


57 802 




88 955 


53 758 


91641 


18 867 




33181 


72 664 


53 807 


00 607 




67 248 


09 082 


12 311 


90 316 




27 589 


99 528 


14480 


50 961 




79130 


90410 


45 420 


77 757 




25 731 


37 525 


16 287 


66181 




45 904 


75 601 


70492 


10274 




19 976 


04925 


07 824 


76 044 




15 218 


49 286 


89 571 


42 903 




1.0423 


—1.8149 


1.1803 




0.0033 


1.8818 


0.7390 


-0.2736 




1.0828 


—1.1147 


0.2776 


0.1012 




-1.3566 


0.6379 


—0.4428 


—2.3006 




-0.6446 


1.4664 


1.6852 


—0.9690 




-0.0831 


—0.2676 


—1.2496 


—1.2125 




1.3846 


—0.6022 


0.0093 


0.2119 




-1.4647 


—0.0572 


—0.5061 


—0.1557 




-1.2384 


1.4943 


—0.4406 


—0.2033 




-0.1316 


—0.8513 


1.1054 


1.2237 




-0.7003 


0.7165 


0.8563 


—1.1630 




1.8800 


distribution (22 


) (see 3.1) 








variable %, with 


a = 0, a=l 
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DIVIDING A SEGMENT IN A GIVEN RATIO 
N M. Beskin 

One of the series "Popular Lectures in Mathematics", widely 
used in the extracurricular maths clubs and circles of Soviet 
schools. Is broadly accessible to pupils doing "0" level ma- 
thematics. Describes various theories arrived at through 
deeper study of exercises in dividing a segment in a 
given ratio. 

Through looking into this problem and related questions, the 
reader makes a short journey through mathematics, encoun- 
tering on the way affine and projective geometry and group 
theory (generally without their labels). 

THE METHOD OF MATHEMATICAL INDUCTION 
I. S. Sominsky 

One of the series "Popular Lectures in Mathematics", widely 
used in school extracurricular maths circles. Will interest 
sixth-formers doing "A" level mathematics, and first-year stu- 
dents whose courses include mathematics. Describes the me- 
thod of mathematical induction so widely used in various fields 
of mathematics from elementary school courses to the most 
complicated investigations. Apart from being indispensible 
for study of mathematics, the ideas of induction have gen- 
eral educational value, and will interest a broad readership 
not familiar with mathematics. 



